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ABSTRACT 

O . 

O I We discuss the current status of a computational approach which allows to evaluate the 
■ dielectric matrix, and hence electronic excitations like optical properties, including local 
field and excitonic effects. We introduce a recent numerical development which greatly 
reduces the use of memory in such type of calculations, and hence eliminates one of the 
bottlenecks for the application to complex systems. We present recent applications of the 
method, focusing our interest on insulating oxides. 
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INTRODUCTION 



Density-functional theory (DFT) in the Local Density Approximation (LDA) is widely and 
(3 I successfully used as a state-of-the-art tool to compute the ground-state electronic properties 
of many-electron systems |l|] . It is feasible to apply DFT to systems as complex as surfaces, 
I defects, and clusters. However, while the ground-state properties can in principle be obtained 
exactly within DFT, many spectroscopic properties are in general not directly accessible in 
. _ such a calculation. 
! In fact, it is well known that the use of the DFT-LDA eigenvalues as the physical energies 

entering in the absorption process relies on several approximations. First of all, Kohn-Sham 
Q i LDA eigenvalues are only a rough approximation of the electron addition or removal energies, 
such those measured in photoemission. The correct energies should be obtained by using the 
true electron self-energy operator S, which appears at the place of the DFT-LDA exchange- 
^ ■ correlation potential in an equation similar to the Kohn-Sham one. A good approximation 

^ I for E, which allows to determine the quasiparticle (QP) bandstructure, can be obtained 
within Hedin's GW approach 0. In this scheme, DFT-LDA results can be used as the 
starting point, in a first-order perturbative approach 0, ^. 

However, the true QP energies, together with the LDA wavef unctions, are in principle 
still not sufficient to describe correctly an absorption process, in which electron-hole pairs 
are created. Their interaction can lead to bound exciton states which occur within the gap, 
and can also induce appreciable distortions of the spectral lineshape above the continuous- 
absorption edge. 

Since recently |^, H, |, |10|, 0, these excitonic effects can be treated in the ab initio 
framework. Various applications have demonstrated the success of the Green's functions 
approach ||^, |^, ^ ^ [TO, |TT], T^, |TB[. However, the calculations are still very cumbersome, 
since a two-particle problem has to be solved. This prevents the method from being applied 
on a large scale. 

In this work, we introduce a perturbative approach to the solution of the Bethe-Salpeter 
equation which greatly reduces the use of memory in such type of calculations, and hence 
eliminates one of the bottlenecks for the application to complex systems. We show that the 
approximation is perfectly controllable, hence the approach is still ab initio. In order to test 
the limits of the method, we present recent applications for the case of insulating oxides. 



THE APPROACH 



The absorption spectrum is given by the imaginary part of the macroscopic dielectric function 

eM{uj) = 1 - hni t;(q)xG=o,G'=o(q; ^^), (1) 

q^O 

where x{r,r';u) = —iS{r,r,r',r';u). The four-point function S obeys the Bethe-Salpeter 
equation, 

^(1, 1'; 2, 2') = Soil, 1'; 2, 2') + ^o(l, 1'; 3, 3')H(3, 3'; 4, 4')^(4, 4'; 2, 2'). (2) 

The notation (1,2) stands for two pairs of space and time coordinates, (ri, ti; r2, ^2)- Re- 
peated arguments are integrated over. The term 5*0(1, l';2,2') = G(1',2')G(2, 1) yields the 
polarization function of independent quasiparticles XO) from which the standard RPA 
is obtained (G(l, 1') is the one-particle Green's function The kernel H contains two 

contributions: 

S(l, 1', 2, 2') = -^5(1, l')5(2, 2')vil, 2) + z5(l, 2)5(1', 2')iy(l, 1')- (3) 



Considering the first term in the calculation of S is equivalent to the inclusion of local 
field effects in the matrix inversion of a standard RPA calculation. In order to obtain 
the macroscopic dielectric constant, the bare Coulomb interaction v contained in this term 
must, however, be used without the long range term of vanishing wave vector [|14|. When 



spin is not explicitly treated, v gets a factor of two for singlet excitons. In the second 
term, W is the screened Coulomb attraction between electron and hole. It is obtained as a 
functional derivative of the self-energy in the GW approximation, neglecting a term 
This latter term contains information about the change in screening due to the excitation, 
and is expected to be small [|I5|. We limit ourselves to static screening, since dynamical 
effects in the electron-hole screening and in the one particle Green's function tend to cancel 
each other |]l6l, which suggests to neglect both of them. The set of equations (1) - (3) are 
at the basis of all the ab initio exciton calculations which have appeared in the literature 
recently [|, 0, |, |, 0, 0, |12|, ig. 



In order to solve Eq. (2), we rewrite it as an effective eigenvalue problem. 



(n3,n4) 

with 

TT{ni,n2),{n3,ni) ^ f _ TP \X X _ j ( f _ f ] 

■'■■'■exc \^'-'n2 -t^rti ;"ni,n3"n2,n4 ''\Jn2 Jrii) ^ 

j ■ j dri dr[ dra dr'^ ^nA^i) Cal^'i) ^(ri, r'^, rs, r'^) i^n^iri) V'n4(r2)- (5) 



The '?/'n(i") are LDA Bloch functions, with n denoting a band index and a Bloch vector 
k. For the calculation of absorption spectra, we can limit ourselves to transitions with 
positive frequency, i.e. (ni,n2) and (n3,n4) are pairs of one valence and one conduction 
band, respectively (in other words, we consider only the resonant part). Moreover, we build 
up the spectra of optical properties by considering only negligible momentum transfer, hence 
the same k for the valence and the conduction state. Equation (1) reads then: 

eM{uj) = l + \imv{q)y2^ — ^ . (6) 



q^o '^'^ {Ex -co) 




Figure 1: Absorption spectrum of bulk silicon. Dots: Experiment |M. Short-dashed curve: 
calculation, no electron-hole interaction included. Continuous curveTmll exciton calculation. 
Long-dashed curve: zero-order result. Dot-dashed curve: first order result. Vertical bars: 
block cuts. 



The problem 

The advantage of this approach with respect to a straightforward inversion of equation (2) 
is that we have to diagonalize the Hamiltonian (5) only once, and can then construct the 

spectrum for all frequencies. Moreover, the knowledge of the coefficients A^^''^'^^ allows for a 
detailed analysis of the results. However, the matrix to be diagonalized can be very large, 
since the basis set is built up of pairs of states. Moreover, the number of k-points in the 
Brillouin zone sampling needed for the calculation of optical properties is typically at least 
one order of magnitude bigger than that employed in a ground state calculation. Even 
though this number enters the basis linearly, and not quadratically as the number of bands, 
the basis is very large, and the calculations are prohibitive for an application to an extended 
range of energies, i.e. including states far from the Fermi level, or to systems with many 
atoms. 

In order to find a way for overcoming this difficulty, it is worthwhile to analyse how 
excitonic effects alter the absorption spectra. 

We ffrst look at bulk Silicon. The importance of excitonic effects in the absorption spec- 
trum of Silicon have already been shown by Hanke and Sham in a semi-empirical calculation 
0. Recently, ab initio calculations p, ^ have confirmed their findings. For this work, we 
have repeated the calculations of Ref. 0, but using a k-point set for the Brillouin zone 
(BZ) sampling which is shifted off the high symmetry directions, whereas in Ref. special 
points [1^ were used. It has in fact been discussed [jl8[ that the shifted grid gives improved 
results, still, the number of k-points which must be used in order to get reasonable results 
is big; here we have chosen 256 points in the BZ, which, together with the 4 valence bands 
and the 3 lowest conduction bands, leads to a matrix size of 3072x3072. (It is obvious that 
increasing the size of the unit cell or the energy range which is considered would let the 
problem explode rapidly.) 
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Figure 2: Contribution of independent-quasiparticle transitions to excitonic transitions at 
given energy (see text), for bulk silicon 



Analysis of the excitonic effect 

The continuous curve in Fig. ^ shows hence the absorption spectrum of bulk silicon which 
we obtain using the parameters described above. It is in good agreement with experimental 
data [0 (dots in Fig. 1), and with the previous ab initio calculations [§, ^. We can 
now analyse the origin of the excitonic effects. First of all, the fact that excitonic effects 
change the peak positions in Silicon can be misleading: in fact, the density of transition 
energies does virtually not change between the independent quasiparticle and the exciton 
result. This has been found not only for silicon, but seems to hold quite generally for those 
materials where no strongly bound excitons exist (see, for example, [|12], |13|). The effect 
is hence entirely due to a change in oscillator strength, in other words, to the coefficients 

yl^"'^^ from equation (4). A^^'^'' expresses the mixing of different independent quasiparticle 
transitions with energy Ec — into a given transition with energy E\. Fig. 2 shows the 

absolute squared value [A^'^'^-'p for three different A, with Ex chosen above the independent 
quasiparticle gap by 0.5 eV, 5 eV and 10 eV, respectively, in function of the energy of 
the contributing independent-quasiparticle transitions, E^ — E^. One can notice that the 
distributions are extremely sharp; in fact more than 90 % of the weight is contained in an 
energy range around each peak position of less than 0.1 eV, and more than 99 % of the 
weight is contained in a region of less than 1 eV. This suggests at first sight that one should 
divide the spectrum into several parts, and for each energy range carry out the calculations 
using the corresponding limited part of the hamiltonian only. 

The perturbative approach 

In other words, one could hope that a separate diagonalisation of several subblocks, neglect- 
ing the interaction between the blocks, could yield a satisfactory absorption spectrum for all 
energies except those very close to the boarder of the blocks. This is however not true: In 



Fig. 1 we compare the resuh obtained by exact diagonahsation of the full 3072x3072 Hamil- 
tonian (continuous curve) with that obtained by choosing about 40 blocks in a reasonable 
way, i.e. by avoiding to cut at energies close to maxima in the density of states. The result 
is given by the long-dashed curve. The vertical lines indicate the block cuts. The agreement 
of the result obtained in this way with the full calculation is not satisfactory at any energy; 
in particular, when we compare with the result obtained by neglecting the electron-hole in- 
teraction (short-dashed curve) we notice that a big part of the excitonic effect has been lost. 
In order to obtain the full effect, it is hence unavoidable to let the states interact even at 
relatively large energy distances. We can however assume that the interaction decreases with 
increasing energy separation, and try to include the missing contributions in a perturbative 



approach, as also suggested in Ref. pO 



In fact, we can treat that portion of the excitonic hamiltonian, which has been neglected 
up to now by our choice of the blocks, in first order perturbation theory for quasi- continuum 
states. In this way, we obtain the dot-dashed curve in Fig. 1: the agreement with the full 
result is now very good. Note that in this way, we have to diagonalize matrices of a maximum 
size of 100x100, whereas the full calculation requires the diagonalization of a 3072x3072 
matrix. This represents a memory saving of about three orders of magnitude, which will be 
very useful for the study of bigger systems, or in order to be able to improve the k-point 
sampling. It must be stressed that the quality of the resulting spectrum can be checked and 
systematically improved by increasing the block size of the zero-order calculation. 



RESULTS FOR INSULATING OXIDES 



This new approach looks hence extremely promising, and it can be expected that the method 
will work similarly well for systems with the same behaviour as silicon: silicon has in fact 
the peculiarity that screening is strong, hence the electron-hole interaction matrix elements 
W^^'^ are weak. The strong deformation of the spectrum comes from the fact that many 
transitions at energies very close to each other, where 

W^^'^ /[{Ec — Ey) — {Ec' — Eyi)] is big, are interacting. This point is of course very favourable 
for our approach, and it is interesting to explore whether the method can also be applied 
in a less favourable situation. In order to become more severe, we can choose an insulator, 
with low dielectric constant, and in the following we will hence look at two insulating oxides: 
MgO and LiaO. 



MgO 



Excitonic effects in the absorption spectrum of magnesium oxide have been calculated ab 
initio by Shirley et al |T0[, using the Haydock recursion method for a direct solution 
of equation (2). Excitonic effects have turned out to alter the absorption spectrum signif- 
icantly. With a dielectric constant of only about eju = 3, MgO is hence a good candidate 
for our test. We have performed the ground state calculations using pseudopotentials of 
the Trouiller-Martins type |2^, a plane- wave cutoff of 50 Ry, and using 10 special points in 
the irreducible part of the BZ. The theoretical lattice constant is 4.13 A in good agreement 
with the experimental value of 4.21 A. We have then determined the GW corrections to the 
DFT-LDA eigenvalues, using 307 reciprocal lattice vectors and 150 bands. The correction 
to the direct gap turns out to be 2.5 eV, and the corrections are slightly larger at the zone 



boundaries, in good agreement with p3[ and p4|. Our independent quasi-particle absorp- 



tion spectrum, calculated using 256 shifted k-points in the BZ and the first 9 bands, (dotted 
curve in Fig. 3) is close to the result of Ref. [|l^, and in very bad agreement with exper- 
iment (dots in Fig. 3) We have then calculated the full exciton absorption spectrum 
according to equation (6), using the same parameters as for the independent quasiparticle 
spectrum. The result is given by the continuous curve in Fig. 3, and is again close to the 
findings of Ref. [^. As in that work, the very strong excitonic effect significantly improves 
the agreement with experiment |2^, |2H] . The main discrepancy is the overestimation of the 



peak intensities, which might be oue to our limited k-point sampling. Now we want to try 
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Figure 3: Absorption spectrum of magnesium oxide. Dots: Experiment |2^. Dotted curve: 
calculation, no electron-hole interaction included. Continuous curve: full exciton calculation. 
Dashed curve: first order result. Vertical bars: block cuts. 



our perturbative approach. First, we can check if this might be reasonable, by analysing the 
exciton eigenstates as we did in the case of silicon. This analysis is shown in Fig. ^. Also 
in MgO, the distributions are looking rather sharp, but of course, due to the stronger W^^'^ , 
transitions from a wider range of the spectrum must be taken into account. This time, about 
90 % of the weight of the peaks is contained in a range of less than 1.5 eV, and in order to 
retain 99 % of the weight we must go up to almost 7 eV for the second peak (3 eV for the 
first one and 2 eV for the third one). This tells us that we certainly have to choose bigger 
blocks than in Silicon in order to obtain reasonable results. However, still we will try to 
reduce the memory requirement by about two orders of magnitude, and therefore we choose 
the block borders as indicated by the vertical lines in Fig. 5. This corresponds to blocks 
of a maximum size of 500x500. The resulting first-order spectrum is given by the dashed 
curve. Again, the agreement with the full calculation is good, and most of the disagreement 
between the independent quasiparticle and the experimental result has been removed. The 
most important errors occur of course at higher energies, where the states are denser and 
hence the energy range of the blocks smaller. Due to the Kramers Kronig relation, one could 
suspect that this fact might spoil the real part of even at low energies, and hence show up 
in spectra linked to both the real and the imaginary part, like reflectance. We have therefore 
used our results in order to calculate the reflectance spectrum of MgO. The independent 
quasiparticle result is given by the dotted curve in Fig. ^ and the result of the full exciton 
calculation is given by the continuous curve. As in the case of Im{eM), the inclusion of the 
electron- hole interaction has a dramatic effect. In particular, the weak structures up to 15 
eV in the GW spectrum are transformed into pronounced peaks, which are consistent with 
the experimentally found structures in that energy range in the reflectance spectrum (see 
inset) p5[. The first-order result, for the block size of 500x500, is given by the dashed curve. 
Again, the agreement with the full calculation at lower energies is excellent, and also the 
structures at higher energies are considerably improved with respect to the GW calculation. 
The performance of the method for MgO is hence less spectacular than in the case of silicon, 
but still the approach turns out to be valid. 
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Figure 4: Contribution of independent-quasiparticle transitions to excitonic transitions at 
given energy (see text), for magnesium oxide 
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Figure 5: Reflectance spectrum of magnesium oxide. Inset: Experiment Dotted curve: 
calculation, no electron-hole interaction included. Continuous curve: full exciton calculation. 
Dashed curve: first order result. 




Figure 6: Absorption spectrum of lithium oxide. Short-dashed curve: calculation, no 
electron-hole interaction included. Continuous curve: full exciton calculation. Long-dashed 
curve: first order result, block size 300. Dot-dashes curve: first order result, block size 700. 
Vertical bars: block cuts 



U2O 

Another insulator of interest for our purpose is lithium oxide, Li20. Li20 is a material 
which is experimentally studied essentially because of its potential importance in the atomic 
energy domain. Several spectroscopic measurements are published in the literature p^, often 
focusing on the characterisation of radiation defects. Recently, there have been theoretical 
and experimental ||2^ studies devoted to the question of excitonic effects at the band 
edge of Li20. In Ref. |Ja lowering of the optical gap due to excitonic effects of the order 
of 1 eV was predicted. The experiment of Ref. seems to be in agreement with this 



finding. It is hence interesting, on one side, to study excitonic effects in Li20 not only at the 
absorption onset, but over a 20 eV range, and to compare with experiment. On the other 
hand, the fact that now also the transition energies are heavily affected by the electron-hole 
interaction presents an additional test for our approach. We have performed the ground state 
calculations using pseudopotentials of the Trouiller- Mart ins type a plane- wave cutoff of 



80 Ry, and using 2 special points in the irreducible part of the BZ. The theoretical lattice 
constant is 4.534 A in good agreement with the experimental value of 4.573 A. We have then 
determined the GW corrections to the DFT-LDA eigenvalues, using up to 541 reciprocal 
lattice vectors and 280 bands. The correction to the direct LDA gap of 5.6 eV turns out to 
be 2.3 eV, again with corrections which are slightly larger at the zone boundaries, in good 
agreement with Ref. 

Our independent quasiparticle absorption spectrum, calculated using 256 shifted k-points 
in the BZ and the first 9 bands, (short-dashed curve in Fig. consequently shows an onset 



at about 7.9 eV, which well compares with the value of 7.99 eV estimated in Ref. p8[ for 
the band-gap energy, and is considerably higher than the measured onset at about 6.5 - 7 
eV ]2^. We have then calculated the full exciton absorption spectrum according to 



equatiori(6)7using the same parameters as for the independent quasiparticle spectrum. The 
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Figure 7: Reflectance spectrum of litliium oxide. Inset: experiment |^8[. Sliort-daslied 
curve: calculation, no electron-hole interaction included. Continuous curve: full exciton 
calculation. Long-dashed curve: first order result, block size 300. Dot-dashed curve: first 
order result, block size 700. 



result is given by the continuous curve in Fig. The onset has now been lowered, to about 7 
eV, which agrees well with the experiment and with the prediction of Ref . . Moreover, as 
in the case of silicon and MgO, considerable oscillator strength has been transferred to lower 
energies, giving rise to a sharp peak at the onset. This behaviour is similar to what has been 
found for LiF in Ref. [|1^]. Our calculated reflectance spectrum, given by the continuous 
curve in Fig. ^ reproduces the main features of the experimental result (see inset Fig. 

f), with a first strong peak at 7 eV, a second broader peak above 8 eV and a sharp minimum 
etween the two structures. A detailed analysis will be necessary in order to check whether 
the small structure on the low energy shoulder of the 8 eV-peak is related to the shoulder 
found experimentally, and whether the excitonic origin of this structure, proposed in Ref. 
2^, can be confirmed. 

For this work, we concentrate on the question whether our perturbative approach can 
also be applied to this system with large excitonic binding energy. We have therefore also 
calculated the first order spectrum, given by the long-dashed curve in Fig. ^. Here the 
block size has been chosen to be 300x300. The essential part of the excitonic effect is 
already included in this spectrum. In particular, the binding energy of the first peak is given 
correctly. As mentioned above, we can improve this result systematically by increasing the 
block size: The dot-dashed curve has been calculated with a block size of 700x700. The 
vertical bars are the corresponding block cuts. The result is now in very good agreement 
with the full one. The same holds for the reflectance spectrum, as can be seen from the 
short-dashed (GW), long-dashed (300x300) and dot-dashed (700x700) curves in Fig.0. This 
means that, even in a very unfavourable situation, our approach is advantageous. 



OUTLOOK 



We should stress the fact that up to now we have concentrated our effort only at the aim to 
save memory. It is however obvious that we can in principle use the same method also to 



save considerable CPU time and disk storage requirements: the biggest part of the exciton 
Hamiltonian W^^^ is now not longer used in a direct diagonalization, but, for the calculation 
of Ax, in the perturbation formula 

Hr,V{E, - E,)Af\ with Htr = Af)'~^^)Hi:^'^''^Af'^^''\ The fact that now 

H acts on a vector, instead of being inverted or diagonalized, allows to change space for 
each part of the Hamiltonian separately. In particular, one can go back to real space for 
the electron- hole attraction part W , which, according to Eq. (3), is then diagonal in two of 
the four coordinates. In this way, our approach will hence combine both the advantages of 
the effective two-particle picture used in ||^, ^ |11|, and of the recursive inversion approach 
proposed in p^ ]. 

CONCLUSIONS 

In conclusion, we have analysed in detail the nature of the excitonic effects in different 
materials. This has allowed us to design a method which improves the efficiency of the 
ab initio calculation of optical spectra including the electron-hole interaction. By using a 
perturbative approach, the method avoids the diagonalisation of large matrices and therefore 
has the advantage of significant memory saving. Still, it is ab initio in the sense that one 
can converge the results smoothly to the exact ones, by a controllable parameter (i.e. the 
block size). 

The method works extremely well for systems with weak interaction, but strong excitonic 
effects, as we have shown for the case of bulk silicon. We can therefore expect that it will be 
of great usefulness for the study of systems like those, for example, built up of many silicon 
atoms (e.g. amorphous silicon), where we are in the ideal condition of weak interaction but 
large density of transitions. 

We have tested the approach also for unfavourable cases, using as examples two insulators. 
In spite of the strongly increased interaction strength, we can still obtain good spectra in 
the gap region, with memory savings of two orders of magnitude. We have also shown 
that the method yields valid results for the real part, not only the imaginary part, of eM, 
therefore allowing for the description of spectra like the reflectance. The inclusion of first 
order corrections to a given subset of states is hence a very promising method in order to 
obtain, with little effort, improved spectra in the low energy region. 
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